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A wide variety of methods have been used to compute percolation thresholds. In lattice percolation, 
the most powerful of these methods consists of microcanonical simulations using the union-find 
algorithm to efficiently determine the connected clusters, and (in two dimensions) using exact values 
from conformal field theory for the probability, at the phase transition, that various kinds of wrapping 
clusters exist on the torus. We apply this approach to percolation in continuum models, finding 
overlaps between objects with real-valued positions and orientations. In particular, we find new values 
of the percolation transition for disks, squares, rotated squares, and rotated sticks in two dimensions, 
and confirm that these transitions behave as conformal field theory predicts. The running time and 
memory use of our algorithm are essentially linear as a function of the number of objects at criticality. 



What happens to the hole when the cheese is gone? 

Bertold Brecht, Mother Courage 



I. INTRODUCTION 

For more than 50 years, percolation theory has been 
used to model static and dynamic properties of porous 
media and other disordered physical systems Ol-ll]- 
Most natural systems correspond to continuum perco- 
lation, yet most analytical and numerical work has fo- 
cused on lattice percolation. This is reasonable since 
continuum and lattice percolation lie in the same uni- 
versality class. For properties that are non-universal, 
however, such as the location of the threshold, one has 
to study discrete and continuum models individually, 
and it is also satisfying to confirm universality experi- 
mentally by measuring critical exponents and crossing 
probabilities. 

In this contribution we discuss an algorithm to com- 
pute the location of the transition in continuum perco- 
lation models. The algorithm works in arbitrary dimen- 
sions, and for arbitrarily shaped objects; here we focus 
on two-dimensional percolation with disks, squares that 
are aligned or randomly rotated, and randomly rotated 
sticks (see Figured]). Our algorithm is an adaption of the 
union-find algorithm of Newman and Ziff [4], the fastest 
known algorithm for lattice percolation. We show that it 
can be adapted to continuum percolation with the aid of 
some simple additional data structures, and we back up 
our claim by computing numerical values of the transi- 
tion points that extend the accuracy of previously known 
values by several orders of magnitude. 

In two-dimensional continuum percolation, a number 
n of penetrable objects are thrown at random in a square 
of size L^. If the mean density p = n/L^ is finite as n and 
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Figure 1. Continuum percolation with disks, randomly rotated 
sticks, and aligned or rotated squares. In each example, the 
wrapping cluster is marked by color. 



L go to infinity, the spatial distribution of the objects' 
centers is a Poisson point process with density p. The 
system percolates if there exists a cluster of overlapping 
objects that spans the square. We follow [4] in using 
periodic boundary conditions, and focusing on clusters 
that wrap around horizontally, vertically, or both. These 
wrapping clusters display better finite-size effects than 
crossing clusters on open boundary conditions. 

If each object has area a, then the probability that a 
percolating cluster exists in the limit L ^ oo clearly de- 
pends only on the product rj = pa. This dimensionless 
quantity is called the filling factor. It also gives the total 
fraction (p of the plane covered by the objects. 



= 1 - e"^ . 
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While we expect continuum percolation to be in the same 
universality class for any fixed shape, the location of the 
transition, i.e., the critical filling factor r]c, depends on the 
shape of the objects. We write rj^, 77°, rjf, and t]^ for the 
percolation of disks, aligned squares, randomly rotated 
squares, and randomly rotated sticks. In defining 77, we 
treat sticks of length ^ as if they have area a = f-. 

Table I] lists the most accurate numerical values for 
Tic from previous work and the work presented here. 
The best previous results on disk percolation are due 
to Quintanilla, Torquato, and Ziff 13], who varied the 
density of the Poisson process as a function of position 
and kept track of the front of the connected cluster. The 
best previous results on squares are due to Baker et al. 10. 
The best previous results on sticks are due to Li and 
Zhang 10], who used an approach similar to ours but 
with open boundary conditions. 

Our results are consistent with the rigorous bounds 

1.127 < r]? < 1.12875 

I. 098 < r]° < 1.0995, (2) 

computed with 99.99% confidence by Balister, Bollobas 
and Walters [8] using a Monte Carlo estimate of a high- 
dimensional integral. On the other hand, it is a little sad 
to dash the hope — which one might have entertained 
after reading |6], and which is just barely consistent 
with Q— that (/)° is exactly 2/3. 

In the following sections, we review the union-find 
algorithm of 0], how it finds wrapping clusters in peri- 
odic boundary conditions, and how we extend it to the 
continuous case. We show that the running time of our 
algorithm is essentially linear in the number of objects, 
i.e., linear in L^. In addition to estimating the threshold, 
we also measure the finite-size exponent v, giving strong 
evidence that these continuum models are in the same 
universality class as lattice percolation. Finally, we find 
that the probability of a wrapping cluster at criticality is 
precisely that predicted by conf ormal field theory. 

II. THE ALGORITHM 

We will simulate percolation in the microcanonical en- 
semble, i.e., where the number n of objects in the square is 
fixed. In each trial, we add one object at a time, stopping 
as soon as a percolating cluster appears. Following ft 
we keep track of the connected components at each step 
using the union-find data structure. In union-find, each 
cluster is represented uniquely by one of its members. 
We have access to two functions: f ind(z), which finds the 
representative r(z) of the cluster to which object i belongs, 
and merge(z,;), which merges z's cluster and j's cluster 
together into a single one with the same representative. 

Internally, union-find works in a very simple way. 
Each object z is linked to a unique "parent'' ip{i) in the 
same cluster, except for the representative which has no 
parent. When we call f ind(z), it follows the links from z 
to its parent ^(z), its grandparent p(p(z)), and so on, until 




Figure 2. When we call f ind(z), we split and shorten the path 
from i to its representative r{i) by setting the parent of each 
object along the path to be its grandparent. This turns a path 
of length I into two paths of length 111. 



it reaches z's representative r(z). Similarly, merge(z, ;) uses 
f ind(z) and f ind(;) to obtain r(z) and r(;), and declares 
one of them to be the parent of the other, unless r(z) = r(;) 
and they are already in the same cluster. 

The running time of find(z) is proportional to the 
length of the path from z to r(z). If merge(z,;) sensi- 
bly links the smaller cluster to the larger one, setting 
p(r(z)) = r(;) whenever z's cluster is smaller than j's, a 
simple inductive argument shows that these paths never 
exceed log2 n in length. However, we can make these 
paths even shorter using a trick called path compression. 
Since r(z) is the representative of every object ; along the 
path from z to r(z), we can set = r(z) for all of them, 
linking them directly to their representative so that find 
will work in a single step the next time we call it. 

As a result, the amortized cost of the find and merge 
operations — that is, the average cost per operation over 
the course of many operations — is nearly constant. 
Specifically, it is proportional to a(n), when a is the in- 
verse of the Ackermann function 13]. The Ackermann 
function grows faster than any primitive recursive func- 
tion, i.e., any function that can be computed with a fixed 
number of for-loops: faster than an exponential, an iter- 
ated tower of exponentials, and so on lHol] . As a conse- 
quence, a{n) grows incredibly slowly, and the smallest 
value of n such that a{n) > 4 is so large that it can only 
be written with exotic notation. Thus the total running 
time for n objects is essentially 0(n). 

In our implementation, we employ a form of path 
compression that is faster and almost as effective: we 
link each object ; on the path to its grandparent, setting 
= Pipij))- This is known as path splitting, since it 
turns a path of length £ into two paths of length ^/2, or 

+ l)/2 and - l)/2 if f is odd, as shown in Figure 121 
It has the advantage of requiring only one pass along 
the path, and it takes just one line of code. Like path 
compression, it guarantees an amortized running time 
of 0{a(n)) [11]. 

For lattice percolation as in [4], each time we add a 
new occupied site, we can check which of its neighbors 
are occupied, and merge them together with the new 
site. In the continuous case, we have more work to do: if 
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1.128085(2) 
1.12808737(6 


>) 


1.098(1) 
1.09884280(9 




0.9819(6) 
0.9822723(2) 


5.63726(2) 
5.6372858(6) 


previous 
our work 


0.6763475(6) 
0.67634831(2) 


0.6666(4) 
0.66674349(3) 


0.6254(2) 
0.62554075(8) 


0.9964373(1) 
0.996437475(2) 



Table I. Numerical values of critical filling factors 7]^ and area factors c^^ = 1 ~ e '^^ in continuum percolation for disks, aligned 
squares, randomly rotated squares, and randomly rotated sticks. Previous estimates are from (1-0] • 




Figure 3. We divide the plane into square bins whose width 
equals the diameter of the disks. Each disk in a given bin 
(dashed) can only intersect with other disks in the same bin, or 
in the eight neighboring bins. 



we add a new disk (say), we have to find which nearby 
disks it intersects. To do this efficiently, we divide the 
plane into square bins as shown in Figure |3l Each disk 
belongs to whichever bin its center lies in. The width of 
each bin is the diameter of the disks, so that a disk in a 
given bin can only intersect with other disks in that bin 
or the eight bins in its neighborhood. 

On average, the number of disks in each bin is a con- 
stant proportional to p, so we can find all the disks inter- 
secting with each new one in constant time. We use the 
same approach for the other shapes; for rotated squares 
of width ^, the bins need to have width ^/2£. A similar 
approach for rotated sticks was used in 

If we wished to detect crossing clusters — those that 
connect, say, the top and bottom edges of the square — 
we could add two special objects to the union-find data 
structure, which are connected by fiat to all the disks in 
the bins along the top or bottom edge. We would then 
check, at each step, whether these two objects are in the 
same cluster. However, as discussed below and in 10], 
the finite-size scaling is much better if we use periodic 
boundary conditions instead, and look for clusters that 
wrap around the torus horizontally or vertically. 

We detect these wrapping clusters using a technique 
originally used for detecting crossing clusters in the Potts 



model (T3]. We associate a vector with each object in the 
union-find data structure, recording the displacement 
between it and its parent. In principle this displacement 
is real-valued, but it suffices to record an integer vector 
giving the displacement between their respective bins. 
When we compress and splint a path, we sum these 
vectors to get the total displacement between each object 
on the path and its new parent. 

Now suppose that merge(z,;) finds that two overlap- 
ping disks i and ; are already in the same cluster. Object 
i now has two paths to its representative; one that goes 
through its own parent, and another that consists of hop- 
ping to ; and then going through fs parents. We sum 
the displacement vectors along both these paths. If these 
sums are the same, then the cluster is simply-connected. 
But if they differ by +L in either coordinate, then the 
cluster has a nontrivial winding number around one or 
both directions on the torus. 

Like the union-find algorithm itself, the time it takes 
to sum these vectors is proportional to the length of the 
paths from i and ; to tlieir representative. As Figure IH 
shows, the total running time of our entire algorithm — 
the time it takes to carry out a trial on a lattice of size 
L, adding objects one at a time until a wrapping cluster 
appears — is essentially linear in the number n of objects 
at criticality, or equivalently linear in L^. It slows down 
somewhat when the computer is forced to switch to parts 
of its memory with slower access, but this only affects 
the leading constant. 



III. ANALYSIS AND RESULTS 

If in each trial we stop at the first n where a wrapping 
cluster appears, then the estimated probability Fi{a,n) 
that a wrapping cluster exists in the microcanonical en- 
semble with n objects of area a is the fraction of trials that 
stop on or before the nth step. To obtain the probabil- 
ity Rhi^) of percolation in the grand canonical ensemble 
with filling fraction 7], we convolve Pi with the Poisson 
distribution with mean A = pV- = rjL^/a: 



RLin) = e-^]_^ — PL{a,n). 



(3) 



n=0 



To avoid numerical difficulties where the numerator and 
denominator are both very large, we compute Poisson 
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Figure 4. Average CPU time T for a single realization of an 
L X L lattice up to the percolation transition. T is measured on 
a laptop with Intel Core 2 Duo 2.53 GHz CPU with 3 MB cache. 
The dashed lines T oc are guides for the eye. The slope 
increases when the cache memory is exhausted, forcing the 
computer to switch to regions of memory with slower access; 
however, the running time remains linear inn ocL^. 



w^eights Wn oc A"/n! inductively in two sequences Wn-k 
and Wfi+k to the left and right of the peak at n = [A J, 
where we define Wn = 1: 



Wn-k = 



I n-{k-l) 



W 



n-{k-l) 



f or fc = 
for fc = 1,2,. 



and 



f or fc = 



^"""^ 1 iTfc '^fi+k-i f or fc = 1, 2, . . . 

The sum (|3|) only needs to be computed for a finite num- 
ber of terms. In one direction, we only need to sum down 
to the smallest n where Pl(^, n) is nonzero, i.e., the small- 
est value of n where we observed a wrapping cluster in 
at least one trial. In the other direction, once we pass 
the largest n where a wrapping cluster first appeared, 
then Pi{a,n) = 1. At that point, we sum the remaining 
terms until they are zero to within the numerical preci- 
sion of the computer. We then normalize the entire sum 
by dividing by X, ^n- 

Equipped with the data from the microcanonical sim- 
ulations and this convolution routine, we compute the 
wrapping probability functions for various system 
sizes L and shapes. Like (3], we look for several kinds of 
wrapping in particular. Specifically: 

• R\(vi) is the probability of any kind of wrapping 
cluster. This is indicated by a winding number 
that is nonzero in either coordinate. 



Figure 5. Wrapping probabilities R\(ri) for disk percolation and 
L = 16, 32, 64, 128, 256, 512. The dashed line is the exact value of 
the critical wrapping probability Rloi^c) from conformal field 
theory. 



R^{7]) is the probability of a cluster that wraps hor- 
izontally. This is indicated by a winding number 
that is nonzero in the first coordinate. 



• Rl(r]) is the probability of a cluster that wraps both 
horizontally and vertically. This is indicated by a 
single winding number that is nonzero in both co- 
ordinates, or a pair of winding numbers that are 
nonzero in the first and second coordinates respec- 
tively. 

• R[{r]) is the probability of a cluster that wraps hor- 
izontally, both not vertically. This is indicated by 
a winding number that is nonzero in only the first 
coordinate. 

For any L and any 77, these probabilities obey 
Rl = 2Rl-Rl = 2Rl + Rl 

We assume here that the torus is square, so that horizon- 
tal and vertical wrapping probabilities are equal. 

Note that if the first nonzero winding number ob- 
served in a given trial is nonzero in both coordinates, 
then a cluster of type 1 (horizontal but not vertical) does 
not occur at all in that trial. Thus R^ir]) does not tend to 
1 as 77 increases. 

In practice, we focused on and R^^. In each run, we 
recorded the number of objects at which horizontal 
wrapping first occurred, and the number at which 
vertical wrapping first occurred. Then = min(n^, n^) 
and = max(n^, n^) are our estimates, in that run, of the 
n at which i^^ and i^^ jump from to 1. 
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Figure 6. Slope of R^iri) at r]L, the estimated critical filling factor. 
The line is 0.361L^/^, confirming the universal critical exponent 
V = 4/3 for finite-size scaling. 



Figure 7. Convergence with increasing system size L of Rl{t]^) 
to its known value at L = oo for disk percolation. The line is 
proportional to L~^, the conjectured convergence of Rl(ric). 



A beautiful fact is that, even though the percolation 
threshold rjc is not know^n for any of our models, conf or- 
mal field theory implies exact values for these probabil- 
ities at the transition in the limit L ^ oo 13, |l^. Specifi- 
cally, 

Rt = 0.521 058 289 248 821 787 848... 
Rl, = 0.690473 724570168677230... 

(4) 

Ri = 0.351642853 927474898465... ^ 
Ri = 0.169415435321346889383... 

For each L, and each type of wrapping cluster, we can 
estimate the critical filling factor rji as the solution of the 
equation 

RM = Roo • (5) 

For instance, Figure|5]shovv^s R^it]) for disks for L ranging 
up to 512. The filling factors rji where these curves cross 
Rl^ rapidly converge to rjc. 

The rate of convergence is determined by two factors. 
The first comes from the fact that the width of the tran- 
sition window from ~ to I^l ~ 1 scales as L~^/^ 
where v = 4/3 is a universal critical exponent for two- 
dimensional percolation. This scaling holds even for 
small systems, as can be seen in Figure [6l where we plot 
the slope of Ri at the estimated critical filling factor rji. 
The slope scales perfectly like L^^^. 

The second factor comes from the fact that Rl(^) not 
only becomes steeper but also moves upward in the crit- 
ical region (see the inset in Figure O. To measure the 
contribution from this effect, we computed the difference 
Rii^c) - Rio using the previously best known value for 
rjc from Table [D This difference scales like L~^, as can be 
seen from Figure [71 The same scaling has been observed 




Figure 8. Estimated critical filling factors for continuum per- 
colation of squares, derived from (jSj with (top) and 
(bottom). Note the resolution of the r]-axis. 

in two-dimensional lattice percolation 0]. Note that the 
periodic boundary conditions are responsible for this 
quadratic decay. With open boundary conditions, this 
factor scales as |T^, leading to more severe finite-size 
effects. 

These two factors combine to give 

nL-ilc-L-^-'"' = L-''" (6) 

for the rate of convergence. Hence we expect a straight 
line if we plot rji vs. L~^^/^, and this is exactly what we 
observe in Figure [HI Extrapolating this line to zero then 
gives our estimates of rjc shown in Table [J 
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How do we compute the error bars in our estimates of 
rjc? First consider the fluctuations in Riit]). Each of our 
microcanonical experiments contributes to our estimate 
of Rl{t]) for all rj through the convolution (|3|. We can 
imagine this as choosing n from the Poisson distribution, 
adding n objects, and returning an estimate of Riir]) = 1 
or depending on whether they percolate or not. If 
we perform N trials, the number of trials that return 1 is 
binomially distributed with mean Ri (r])N, and averaging 
gives an estimate of Riit]) with standard deviation 

Depending on which kind of wrapping cluster we are 
looking for, this is roughly OAN~^^^. 

When we look for the rji where Rl{7]) crosses Roor the 
error on rji is given by 

Since the slope R[{riL) grows as 0.361L^/^ (see Figure O 
this gives 

These are the error bars shown in Figure O 

The extrapolated value for rjc is computed from simu- 
lations for D different system sizes L, which in a weighted 
linear regression as in Figure [8] yields an error roughly 
^fD times smaller than the error bars of the underlying 
data points. 



Finally, we average our estimates of rjc from R^^ and 
i^^. Assuming that these estimates are only weakly cor- 
related reduces the error bars by another factor of V2. 

The error bars shown in Table [J are the result of sim- 
ulating roughly D = 50 system sizes ranging from L = S 
to L = 2048, with sample sizes N ranging from 10^° for 
the systems with L < 100, to 10^ for 100 < L < 500, to 10^ 
for 500 < L < 2048. 

We ran these simulations in parallel on several 
computer clusters with greatly varying computational 
power. In total, our simulations would have taken about 
400 years if done only on the laptop on which this paper 
was written. 



IV. CONCLUSIONS 

We have shown that the union-find approach to esti- 
mating percolation thresholds introduced by Newman 
and Ziff [4] can be applied in the continuous case. With 
the help of an algorithm for estimating rjc that runs in 
essentially linear time as a function of the number of ob- 
jects at criticality, we have obtained new estimates for 
rjc in a variety of continuum percolation models that are 
several orders of magnitude more accurate than previous 
results. In the process, we have confirmed the predic- 
tions of conformal field theory for these models, both for 
the finite-size scaling exponent v and the probabilities 
that various kinds of wrapping clusters exist at rjc on 
periodic boundary conditions. 
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